%let mNumScenarios = 10; 
/* number of scenarios per period */

proc optmodel;
	num nIns = 1;
	/* nIns is the number of instances */
	set<num> INSTANCES = 1..nIns;
	num seed = 1131;

	num price_lb{INSTANCES} = .5;
	num price_ub{i in INSTANCES} = 1.5;
	num salvage{i in INSTANCES} = .3;
	num nPeriods{INSTANCES} = 3;
	num seed_per_ins{INSTANCES} = round(ranuni(seed) * 10000) - 1;
	/* nPeriods is the number of periods */

	create data main from [instance]=INSTANCES nPeriods price_lb price_ub salvage seed=seed_per_ins;
quit;

%macro create_data(nPeriods);
	data scenario_tree_size;
		do period=1 to &nPeriods;
			nsamples=&mNumScenarios;
/*			nsamples=round(exp(log(&mNumScenarios)/(&nPeriods-1)));*/
			output;
		end;
	run;

	%let param_factor = ;
	%do i=1 %to &nPeriods;
		data _null_;
			set scenario_tree_size(firstobs=&i obs=&i);
			call symputx('param_factor',"&param_factor "||compress('period_'||period||'='||nsamples)||' ordered');
		run;
		%put &param_factor;
	%end;
	proc plan;
		factors &param_factor;
		output out=sample_path;
	quit;

	proc optmodel;
		set PERIODS = 1..&nPeriods;

/*receipts*/
		num receipt{PERIODS} init round(3*ranuni(&seed));
		receipt[1]=12;
		receipt[2]=12;
		receipt[3]=12;
		create data receipts from [period]=PERIODS receipt;

		set<num> SCENARIOS;
		num path{SCENARIOS, PERIODS};
		read data sample_path into SCENARIOS=[_n_] {t in setof{i in PERIODS} ('period_'||i)}<path[_n_,input(scan(t,2,'_'),best8.)]=col(t)>;

		str decision{s in SCENARIOS, t in PERIODS} init '';
		set<str> DECISIONS{t in PERIODS} init {};
		for {s in SCENARIOS, t in PERIODS} do;
			for {i in PERIODS: i<=t} do;
				decision[s,t] = decision[s,t]||'_'||path[s,i];
			end;
			DECISIONS[t] = DECISIONS[t] union {decision[s,t]};
		end;
	
/*demand function*/
		num var_eta = 1;
		num var_epsilon = 1;

		num eta{t in PERIODS, d in DECISIONS[t]}	= exp(-log(var_eta+1)/2 + sqrt(log(var_eta+1))*rannor(&seed));
		num epsilon{t in PERIODS, d in DECISIONS[t]} = 0+sqrt(var_epsilon)*rannor(&seed);
		num eta_per_scenario{s in SCENARIOS, t in PERIODS} = round(eta[t,decision[s,t]],0.01);
		num epsilon_per_scenario{s in SCENARIOS, t in PERIODS} 
			= round(max(-eta_per_scenario[s,t]*(16-4*&price_ub),epsilon[t,decision[s,t]]),0.01);
		create data eta from [scenario]=SCENARIOS {t in PERIODS}<col('period_'||t)=eta_per_scenario[scenario,t]>;
		create data epsilon from [scenario]=SCENARIOS {t in PERIODS}<col('period_'||t)=epsilon_per_scenario[scenario,t]>;

		create data postponed_price_decision from [scenario]=SCENARIOS {t in PERIODS}<col('period_'||t)=decision[scenario,t]>;

		for {s in SCENARIOS, t in PERIODS} do;
			decision[s,t] = '_'||s||'_t'||t;
		end;
		create data det_decision from [scenario]=SCENARIOS {t in PERIODS}<col('period_'||t)=decision[scenario,t]>;

		for {s in SCENARIOS, t in PERIODS}	decision[s,t]='';
		for {s in SCENARIOS, t in PERIODS: t < card(PERIODS)} do;
			for {i in PERIODS: i<=t} do;
				decision[s,t+1] = decision[s,t]||'_'||path[s,i];
			end;
		end;
		for {s in SCENARIOS} decision[s,1] = '_';
		create data price_decision from [scenario]=SCENARIOS {t in PERIODS}<col('period_'||t)=decision[scenario,t]>;
	quit;

	proc sort data=sample_path; by 
	%do i=1 %to &nPeriods;
		period_&i
	%end;
	;
	quit;

	proc sql;
		create table solution_postponed_profit
		(
			id 			num 	label='Instance ID',
			profit_ub	num 	label='Profit upper bound',
			profit		num 	label='Profit',
			profit_lb	num 	label='Profit lower bound',
			gap			num 	label='Gap(%)'
		);
		create table solution_delayed_profit
		(
			id 			num 	label='Instance ID',
			profit_ub	num 	label='Profit upper bound',
			profit		num 	label='Profit',
			profit_lb	num 	label='Profit lower bound',
			gap			num 	label='Gap(%)'
		);
		create table data_size
		(
			id 			num 	label='Instance ID',
			type		char(16) label='problem type',
			period		num		label='# periods',
			Branch		num		label='# branch',
			Scenario	num		label='# scenario',
			ncol		num		label='# col',
			nrow		num		label='# row'
		);
	quit;
%mend create_data;

%macro solve_type(problem_type);
	put "solving &problem_type pricing problem...";
	problem_type = "&problem_type";
	put problem_type=;

	if problem_type="POSTPONED" then restore profitlb1;
	if problem_type="DET" then restore profitlb2;

	primalinModel=if count("&problem_type",'_') then scan("&problem_type",2,'_') else '';
	put primalinModel=;

	for {s in SCENARIOS, t in PERIODS: primalinModel ne ''} 
		feaPrice[s,t,"&problem_type"] = 1/card(SCENARIOS_WITH_DECISION["&problem_type",decision[s,t,"&problem_type"]]) *
				sum{ss in SCENARIOS_WITH_DECISION["&problem_type",decision[s,t,"&problem_type"]]} optPrice[ss,t,primalinModel];

	for {s in SCENARIOS, t in PERIODS: primalinModel ne ''} fix Price[decision[s,t,"&problem_type"]] = feaPrice[s,t,"&problem_type"];

	solve obj profit with nlp/ printfreq=20;* /multistart seed=42 msmaxstarts=3;* /tech=ipfilter;
	solution_profit["&problem_type"] = profit.sol;


	for {s in SCENARIOS, t in PERIODS}
		optPrice[s,t,"&problem_type"] = Price_s[s,t].sol;

	if primalinModel='' then do;
		create data opt_price_&problem_type from [scenario]=SCENARIOS {t in PERIODS}<col('period_'||t)=price_s[scenario,t].sol>;

		create data data_size_&problem_type from id=&ins type="&problem_type" period=(card(PERIODS)) Branch=&mNumScenarios 
				Scenario=(card(SCENARIOS))
				ncol=(card(DECISIONS["&problem_type"])+2*card(SCENARIOS)*card(PERIODS))
				nrow=(2*card(SCENARIOS)*card(PERIODS));
	end;
	else do;
		create data fea_price_&problem_type from [scenario]=SCENARIOS {t in PERIODS}<col('period_'||t)=price_s[scenario,t].sol>;
	end;

	for {t in PERIODS}
		estPrice[t,"&problem_type"] = sum{s in SCENARIOS} optPrice[s,t,"&problem_type"];

	unfix Price;
	drop profitlb0 profitlb1 profitlb2;
%mend solve_type;

%macro solve_instance(price_lb, price_ub, salvage, seed);
proc optmodel;
	set PERIODS;
	read data scenario_tree_size into PERIODS=[period];
	
	num receipt{PERIODS};
	read data receipts into [period] receipt;

	set<num> SCENARIOS;
	num path{SCENARIOS, PERIODS};
	read data sample_path into SCENARIOS=[_n_] {t in PERIODS}<path[_n_,t]=col('period_'||t)>;

	num eta{s in SCENARIOS, t in PERIODS};
	num epsilon{s in SCENARIOS, t in PERIODS};
	read data eta into [scenario] {t in PERIODS}<eta[scenario,t]=col('period_'||t)>;
	read data epsilon into [scenario] {t in PERIODS}<epsilon[scenario,t]=col('period_'||t)>;

	set<str> TYPES = {'DELAYED','POSTPONED','DET','POSTPONED_DET','DELAYED_DET','DELAYED_POSTPONED'};
	str decision{SCENARIOS, PERIODS, TYPES};
	read data postponed_price_decision into [_n_] {t in PERIODS}<decision[_n_,t,'POSTPONED']=col('period_'||t)>;
	read data price_decision into [_n_] {t in PERIODS}<decision[_n_,t,'DELAYED']=col('period_'||t)>;
	read data det_decision into [_n_] {t in PERIODS}<decision[_n_,t,'DET']=col('period_'||t)>;
	read data postponed_price_decision into [_n_] {t in PERIODS}<decision[_n_,t,'POSTPONED_DET']=col('period_'||t)>;
	read data price_decision into [_n_] {t in PERIODS}<decision[_n_,t,'DELAYED_DET']=col('period_'||t)>;
	read data price_decision into [_n_] {t in PERIODS}<decision[_n_,t,'DELAYED_POSTPONED']=col('period_'||t)>;

	set DECISIONS{type in TYPES} = union{s in SCENARIOS, t in PERIODS} {decision[s,t,type]};

	set SCENARIOS_WITH_DECISION{type in TYPES, d in DECISIONS[type]} = setof{s in SCENARIOS, t in PERIODS: decision[s,t,type]=d} s;
	num feaPrice{s in SCENARIOS, t in PERIODS, type in TYPES} init .;
	num optPrice{s in SCENARIOS, t in PERIODS, type in TYPES} init .;
	num estPrice{t in PERIODS, type in TYPES} init .;
	num showPrice{s in SCENARIOS, t in PERIODS};
	num solution_profit{type in TYPES};

	str problem_type;
	str primalinModel;

	var Inventory{SCENARIOS, PERIODS}  >=0;
	var Sales{DECISIONS[problem_type]} >=0;
	var Price{DECISIONS[problem_type]} >=&price_lb <=&price_ub;
	impvar Price_s{s in SCENARIOS, t in PERIODS} = Price[decision[s,t,problem_type]];
	impvar Sales_s{s in SCENARIOS, t in PERIODS} = Sales[decision[s,t,problem_type]];

	impvar Demand{s in SCENARIOS, t in PERIODS} 
		   = if eta[s,t]*(16-4*&price_lb)+epsilon[s,t]<=0 then 0
			else eta[s,t]*(16-4*Price[decision[s,t,problem_type]])+epsilon[s,t];

	con demand_con{s in SCENARIOS, t in PERIODS}:
		Sales_s[s,t] <= Demand[s,t];

	con lot_sizing{s in SCENARIOS, t in PERIODS}:
		Inventory[s,t] = (if t=1 then 0 else Inventory[s,t-1]) + receipt[t] - Sales_s[s,t];

	max profit = sum{s in SCENARIOS} ( 
					sum{t in PERIODS} (Price[decision[s,t,problem_type]]*Sales_s[s,t]) + &salvage*Inventory[s, max{t in PERIODS} t]
				 ) * 1/card(SCENARIOS);

	con profitlb0: profit>=max(solution_profit['DELAYED_DET'],solution_profit['DELAYED_POSTPONED']);
	con profitlb1: profit>=solution_profit['DELAYED'];
	con profitlb2: profit>=solution_profit['POSTPONED'];
	drop profitlb0 profitlb1 profitlb2;

	%solve_type(DELAYED);
	%solve_type(POSTPONED);
	%solve_type(DET);
	%solve_type(POSTPONED_DET);
	%solve_type(DELAYED_DET);
	%solve_type(DELAYED_POSTPONED);

	restore profitlb0;
	%solve_type(DELAYED);

	num gap{type in TYPES} = if count(type,'_') then 
			round((solution_profit[scan(type,1,'_')] - solution_profit[type]) / solution_profit[scan(type,1,'_')] * 100,.01)	
			else .;

	create data solution_profit from [type]=TYPES solution_profit gap;
quit;

data data_size_det; format type $16.; length type $16; set data_size_det; run;
data data_size_delayed; format type $16.; length type $16; set data_size_delayed; run;
data data_size_postponed; format type $16.; length type $16; set data_size_postponed; run;
proc append base=data_size data=data_size_det; quit;
proc append base=data_size data=data_size_delayed; quit;
proc append base=data_size data=data_size_postponed; quit;
%mend solve_instance;

%macro solve;
proc sql noprint;
	select count(*) into :nIns from main;
quit;

%do ins=1 %to &nIns;
	data _null_;
		set main(firstobs=&ins obs=&ins);
		call symputx('nPeriods',nPeriods);
		call symputx('price_lb',price_lb);
		call symputx('price_ub',price_ub);
		call symputx('salvage',salvage);
		call symputx('seed',seed);
	run;

	%create_data(&nPeriods);
	%solve_instance(&price_lb, &price_ub, &salvage, &seed);

	data solution_profit_&ins;
		set solution_profit;
	run;
	data data_size_&ins;
		set data_size;
	run;
%end;
%mend solve;

%solve;
